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Abstract. We consider heat transport across a harmonic chain connected at its two 
ends to white-noise Langevin reservoirs at different temperatures. In the steady state 
of this system the heat Q flowing from one reservoir into the system in a finite time 
T has a distribution P{Q,t). Wc study the large time form of the corresponding 
moment generating function {c~^'^) ^ g{\) e'^'^^^\ Exact formal expressions, in terms 
of phonon Green's functions , are obtained for both fi{X) and also the lowest order 
correction g{X). We point out that, in general a knowledge of both ^(A) and g{X) is 
required for finding the large deviation function associated with P{Q, r). The function 
fi{X) is known to be the largest eigenvector of an appropriate Fokker-Planck type 
operator and our method also gives the corresponding eigenvector exactly. 
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1. Introduction 

Among the most interesting recent developments in the theory of nonequilibrium 
processes are the so-called fluctuation theorems [H [2|, [3], HI O [6l [7]. These theorems 
make quantitative statements on the probability of negative entropy production in 
nonequilibrium systems. They have been theoretically [HI El [THl [HI [12], [131 [HI [iSl [El [171 
M, [19] and experimentally [20l[2ll[22l[23l[2ll[25l[26l[271[28] studied in a large number 
of systems in various nonequilibrium states. The results have been obtained both in 
the context of transient and steady state phenomena. In the case of nonequilibrium 
steady states in systems carrying heat or particle current, the fluctuation theorems have 
pointed to the importance of the large deviation function (LDF) and the cumulant 
generating function (CGF) [29]. The steady state fluctuation theorem can, in these 
cases, be equivalently stated as a symmetry property of the LDF or the CGF [121 [El [IS] ■ 
Apart from their interest from the point of the fluctuation theorem these two functions 
contain important information on nonequilibrium processes: the LDF gives the precise 
probability of the occurrence of rare events while the CGF contains information on the 
average current in a system as well as all moments [29]. So far there are very few 
examples where either the LDF or the CGF have been exactly computed. The few 
examples include particle transport in exclusion processes [l2l[TH[T5], Brownian motors 
[I8l [T9] , power dissipation and heat transport in single Brownian particles [lEl [301 [SB [32] 
and heat conduction across a quantum harmonic chain [T7] . 

Consider the example of heat conduction through a system coupled to two heat 
baths at different temperatures and Tr and let Q be the heat flowing from the left 
reservoir into the system during a time interval r. Then Q is a stochastic variable with 
a distribution P{Q) and the LDF and CGF for this problem are deflned by the following 
scaling forms, valid for large r: 

/oo 
rfge-^«P(g,r)~(7(A)e^^("). (1) 
-oo 

We will refer to h{q) and fi{X) as the LDF and CGF respectively. For large r the 
term g{X) is a correction to the CGF and can normally be ignored in the saddle-point 
calculation which relates h{q) and yu(A). The saddle-point calculation gives: 

%) = -[MA*) + gA*], ^'{X*) = -q. (2) 

While the LDF and the CGF are normally related by Legendre transformations, there 
are several examples where this relation is known to break down. This happens if the 
function g{X) has singularities in the region of the saddle-point integration. Interestingly, 
in such cases, the CGF might still satisfy the fluctuation symmetry relation, while the 
LDF does not [lEl [321 [33]. Thus we note that if one is interested in the LDF then it is 
in general important to calculate both /i(A) and the leading correction term g{X). Of 
course n{X) is itself of interest since it contains important information on current and 
current-noise properties and relations between response functions. 
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The aim of this paper is to present a formahsm to obtain /i(A) as well as g{X) 
for the problem of heat conduction across a harmonic chain connected to white-noise 
Langevin heat baths. We use the linearity of the problem and show that the problem of 
finding the generating function Z{X) reduces to performing mult i- dimensional Gaussian 
integrations. We are able to find a closed form expression for yu(A), as given by (ITTl) . in 
terms of the phonon-transmission function, a well known quantity in the study of heat 
conduction in harmonic systems. Finding g{X) is more difficult but we are able to also 
express it, as given by ( l59l) . in terms of appropriate phonon Green's functions. For the 
case of a single free Brownian particle we can use our approach to explicitly obtain both 
/i(A) and g{X) and for this case we recover the results of Visco [I6], obtained by solving 
the Fokker-Planck equation. 

The paper is organized as follows. In Sec. ([2]) we define the model that we study, 
make some general remarks on the problem of computing the generating function for 
heat, and briefiy explain our method. In Sec. ([3]) we give the calculation of the CGF 
//(A) while in Sec. (jl]) we give the calculation of the correction term g{X). The example of 
a single Brownian particle, for which explicit results for yu(A) and g{X) can be obtained, 
are considered in Sec. ([5]). Finally we discuss our results in Sec. (El). 

2. Model and general considerations 

We consider a one-dimensional chain of particles with harmonic interactions and 
described by the Hamiltonian: 

N ^ ^ N N 

^ = XI 2"^'^''^ + 2 ^ ^IrnXlXm, (3) 

1=1 1=1 m=l 

where xi, vi and mi are, respectively, the displacement about the equilibrium position, 
velocity and mass of Ith particle and the matrix $ represents the force matrix of the 
system. For the moment we assume that at least one site of the chain is pinned so that 
the centre of mass attains a steady state distribution. However, later we show that the 
results are also valid for a free harmonic chain. The particles 1 and N at the two ends 
— which we refer as left (L) and right (R) respectively — are coupled to white noise 
Langevin heat reservoirs at two different temperatures Tl and Tr respectively. The 
system, described by the variables X"^ = {xi,X2, ■ ■ ■ ,xn) and = {vi,V2, ■ ■ ■ ,viy), 
evolves according to the following equations of motion: 

X = V, MV = -^X --fV + r]{t), (4) 

where M = diag(mi, m2, . . . , mjy) is the mass matrix, the dissipation matrix 7 
has matrix elements 'jij = 5i,j(5j,i7L + Si,NlR) and the noise vector 77 has elements 
rii{t) = Si^irjiit) + 5i,Ar?7/j(t). The variables rji{t),rjR{t) are zero-mean Gaussian white 
noises with correlations given by: 

{Va(t)ria'{t')) = 26a,a' da S(t—t') , where da = jaTa , a,a' = L, R . (5) 

and we have set the Boltzmann constant /c^ = 1. 
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Since the equations of motion are linear and the noise vector rj is Gaussian, 
the probabihty distribution function of the phase space variables U'^ = {X"^ , V^) in the 
nonequilibrium steady state is a Gaussian with mean (U) = and with covariance matrix 
\imt^ao{UU'^). We denote the nonequilibrium steady state distribution by Pss{U). The 
covariance matrix of the ordered harmonic chain was obtained exactly in [3l] . For mass- 
disordered systems the covariance matrix can be expressed in terms of phonon Green's 
functions [351 ES] . The quantity of our interest here is the total amount of heat, Q, 
flowing from one of the reservoirs — say the left (L) — into the system in a given time 
duration r, in the nonequilibrium steady state. This is given by 

Q= f IVLit) - ^LVl{t)Mt) dt, (6) 

Jo 

where Vi{t) evolves according to (j4]), with the initial condition at t = drawn from 
the nonequilibrium steady state distribution. Clearly, Q is a fluctuating quantity 
whose value depends on the initial conditions Uq = U{t = 0) and the noise trajectory 
{ri{t) : < t < r} during any particular realization. Let P{Q, r) denote the probability 
distribution of Q and let Z{X) = {e~^'^) be the corresponding characteristic function, 
where (...) denotes an average over initial configurations as well as over different paths. 

It is useful to consider the restricted characteristic function Z{X,U,t\Uq) = 
{e~^'^)^^ jj where the expectation is taken over all trajectories of the system that evolve 
from a given initial configuration Uq to a given final configuration U in time r. As shown 



m 



Appendix A , this satisfies a Fokker-Planck-type equation: 

drZ{\ U, r\Uo) = CxZ{X, U, r\Uo) , (7) 

with the initial condition Z{X, U, 0\Uq) = 6{U — Uq). The solution of this can formally 
be written down in the eigenbases of the Fokker-Planck operator Cx, and the large r 
behavior is dominated by the term having the largest eigenvalue /i(A), i.e., 

Z{X, U, t\Uo) ~ xiUo, X)^{U, X) exp[r/i(A)] (8) 

where \E'(f/, A) is the eigenfunction corresponding to the largest eigenvalue, i.e., 
Cx'^{U,X) = /i(A)\l/([/, A), and x{Uo, X) is the projection of the initial state onto 
the eigenstate corresponding to the eigenvalue /i(A). We note that for A = 0, 
Z{X = 0, U, t\Uo) is just the phase space distribution at time r. Hence the existence of a 
unique nonequilibrium steady state, which has been proved for this system [36], requires 
that Z{X = 0, f/, r — )• oo|?7o) = Pss{U) and this implies that /i(0) = 0, x{Uo, 0) = 1 and 
^{U, 0) = Pss{U) . Using and the fact that Pss(^) = "^{U, 0), we get for large r: 

Z{X) = JdUoJ dU^{Uo,0)Z{X,U,r\Uo) ~ 9{X) exp[r/i(A)], (9) 

where ^(A) = J dUo *(f/o,0) xiUo,X) J dU ^(f/,A) . 

Note that g{0) = 1. As discussed in the introduction the large deviation function 
h{q = Q/r) = — limT-^oo lnP(Q, t)/t is given by the the Legendre transformation, 

h{q) = -\fi{X*) + X*q\, (10) 
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with X*{q) implicitly given by the saddle point equation /i'(A*) = — g. The above relation 
holds provided that g{X) is analytic along the real A in the region [0,A*], so that g{X) 
can be neglected in the saddle-point calculation as a subleading contribution and the 
contour of integration can be deformed smoothly through the saddle point A*. On the 
other hand, if g{X) possesses any singularity in the region [0, A*], then the contour of the 
integration cannot be deformed smoothly through the saddle point A*, and one needs 
to include the singular part of g{X) in the saddle point calculation. 

The calculation of the LDF thus requires one to compute the CGF /i(A) and the 
leading correction g{X). From the above discussion we see that these can be obtained 
from the largest eigenvalue and eigenvector of an appropriate Fokker-Planck operator for 
this system. This is however very difficult in most cases, including for the model studied 
here. However the linearity of the dynamics and the Gaussian nature of the noise in 
the present problem allow the computation of /i(A) and g{X) using a different approach. 
The basic idea we use is that the variable of interest Q is a quadratic function of the 
initial phase-space configuration Uq and the noise-trajectories {r]{t) : < t < r}, both 
of which are Gaussian distributed variables. Hence the problem of computing (e~^*^) 
reduces to one of doing a multi-variate Gaussian integration. In the following sections 
we present the details. 

We make some remarks on the symmetry property of the CGF. In general, if the 
operator Cx and its adjoint C\ possess the symmetry C\ = Ca^x, then it immediately 
follows that fi{X) = fi{a — A). Even if another operator C'^^ — which is related to Cx by 
a similarity transformation — possesses the symmetry C\ = C'a-x, then also /i(A) has 
the above symmetry. There are some examples of systems with Markovian dynamics 
where the evolution operator satisfies this property but this does not seem to be the 
case for the model discussed here. In fact, even for the simplest case of a single free 
Brownian particle connected to two heat reservoirs [16], the Fokker-Planck operator does 
not possesses the above mentioned symmetry. In this case however, the Fokker-Planck 
operator can be transformed to a Hermitian operator of a quantum harmonic oscillator 
where the potential remains invariant under A — (A/3 — A). We are not aware of such a 
transformation for a system having more than one particle and hence the symmetry of 
/i(A) is a non-trivial one. 

3. Calculation of CGF /^(A) 

Before giving the details of the calculation we first state our main result of this section 
for the CGF which is: 




(11) 



where A/3 = ~ '^l 




(12) 



with G^{u) 



(13) 
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Note that G (u) = G+(— w) = G+*(a;) (* denotes complex conjugate). Using (ITT]) , one 
can immediately verify the expression for the average energy current, 

liin M = _^'(o) = r duTiu). (14) 

which was obtained previously [STJ [361 EH [39]. From ([TT]) . we note that /i(0) = as 
required. It is also evident that ( ITT]) possesses the symmetry /i(A) = /i(A/3 — A). We 
also note that (]TT]) agrees with the classical limit of the result obtained in [T7| for a 
quantum chain. 

We now give the details of the derivation of the expression (1111) . We solve the 
Langevin equations of motion by Fourier transforms. Let us define the finite-time Fourier 
transforms and their inverses as follows: 

{XK), V{iOn),v{uJn)} = - f dt {X{t), V{t),r]{t)}exp{-tujJ) 



T 

oo 



{X{t),V{t),r]{t)}= J2 {X{un),V{un),vM} exp{iunt) 

n=—oo 

with ojn = . (15) 

The Gaussian noise configurations represented by {rj{t) : < t < r}, can now 
equivalently be described in the frequency domain, by the infinite sequence {rjiun) '■ 
n = — oo,...,— l,0,...,oo} of Gaussian random variables having the correlations 

{j]a{uj)riai{oo')) = 26a,a'——— ^[^ + with a,a' = L, R . (16) 

r 

Henceforth for convenience we will drop the subscript n from Un- Taking the Fourier 
transform of (|1]) gives the velocity vector in frequency domain as 

V{uj) = iuG+ ?i + -G+( $AX, - iuMAVr ), (17) 
r 

where AX, = X(r) - X(0), AV; = V{r)- 1/(0) , 
and the Green's function matrix G~^{u) is given by ([T3]) . Since rj ~ 

l/ri/2 and AX,, AK 

have finite variances for large r, it follows that the second term in (fTTI) is 0(l/r^/^) 
smaller than the first. As will be shown in the next section they contribute to the 
correction term g{X). In this section we focus on computing /i(A) and the first term in 
([T7]) then gives: 

In terms of the Fourier transform, the expression in ([6]) for the heat transfer becomes 



oo 



<5(^) = ^ In with Qn = riL{uj)vi{-uj) - -fLVi{uj)vi{-uj). (19) 



n=— oo 



For a chain with at least one pinned site it is easily seen from ( [T8]) that Vi{0) = 
and hence the heat transfer % through the zero-th mode vanishes in ([T9]) . This is 
related to the fact that for a pinned system there is no zero-frequency mode available 
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for transporting energy. Substituting (fTSl) in (HM we write Q(r) = T^^i((?n + Q-n) 
and using the fact rja{—uj) = ?7*(w) we get: 

oo 

Z{X) = (e-^Q) ~ Yl /e^^"«"^"«") , (20) 



where ^„ = (77^(0;), 77^(0;) 



n=l 

T 



where in obtaining the (1, 1)*'^ element above we have made use of the identity 

-lmGl,{u) = uj [7L|G+iMp + 7/?|Gt^MP] , (22) 

which can be proved as follows. From the definition of the Green's function in (fT3l) 
we have [G"*"]"^ — [G~]~^ = 2iuj^. Multiplying both sides by G"*" from the left, and 
by G~ from the right we get G~ — = 2iuG~^ lG~ . The (1, 1)**^ element of this 
leads to ( 1221) . For each n, the average in (l20ll is evaluated with respect to the Gaussian 
distribution 

= ;^mLd (-^" ^''^") ' ^ = ^i^g ( V'' V^) • ^23) 

Hence we get (see [Appendix Bj ) 

e = [det(J + Ari:)A„)J , (24) 

T{uj)TlTr\{^I3-X)\~\ (25) 

where T(cj) is given by f|T2|) . Using (12011251) and noting that in the r — 00 limit we can 
replace the summation over n by an integral over u we obtain our final result in ( ITT]) . 

We remark that, the result (1111) is in fact valid for both pinned as well as unpinned 
cases. For the unpinned harmonic system the centre of mass coordinate does not reach 
a steady state and it is useful to separate out this degree of freedom. Let us define new 
relative coordinates yi = xi — xn, I = 1,2, . . . , N — 1. The facts that the unpinned 
system has translational symmetry and that $ is a symmetric matrix imply the relation 
T.i=i,N = T.j=i,N = 0- Using this we get from dH): 

m =vi-VN, l = l,2,...,N -1 (26) 

Af-l 

mivi = ^ijVj - + Vi{t), l = l,2,...,N . (27) 

i=i 

The relative coordinates yi, I = 1,2, . . . ,N — 1 attain steady state with finite variance 
and it follows then from ( l26l) that viiu = 0) ~ vn{uj = 0) up to order 0(l/r). Then 
from (p7|) we get vi{0) ~ [7/l(0) + //_r(0)]/[7l + 'Jr] for all /. We use this to compute the 
go ^ term in With = (^l(O), ^r(O))^, we get 



% ~ l^^Ao^o where = | %^^t\f ] • (28) 



(7l+7-r)^ (7L+7ii)^ 
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The noise has the Gaussian distribution 

1 

27rVdetD 



and hence we get 



-Argo ^ 



-1/2 



-1/2 



)~ det{I + \tDAo) = 1 + TiO)TLTnX{A(3 - X) (30) 

where T(0) = '^''^lIr/ {il + 7_r)^ is precisely what one obtains by taking T{uj — )■ 0) in 
( !T2|) for the unpinned system. This follows from the fact that in this case G^jy|aj-j.o ~ 
l/iuj{'~)L + 7fl); a result which we now prove: For any matrix A let ^1^'';) denote the 
submatrix of A that occurs between the {i.,jy^ and the (/c, Z)**^ elements. Also let 
B = -Mu;2 + $ + itj^ . Then we have = det Bf^\^j^^ j det Taylor 

(\ 1'] (\ 2) 

expanding the determinants about a; = we obtain det B^j^_^ = det ^) + 0(1^) 

and det -B = 2^7^ det + iaJ7i{ det + O(aj^), where we have used 

det $ = that follows from the property J2i=i n ~ X]j=i n ^iJ ~ 0- Using the latter 
property it is easy to show that det = det jy_-^-) = (— det ^-j. 

Hence we get the desired result. 



4. Calculation of g{\) 

We turn now to the more difficult problem of calculating g{X) which requires one to keep 
the second term in (fT7|l and perform the average over initial conditions, in addition to 
the noise average. We recall that the heat transfer is given by ( lT9l) where the velocity 
Vi can be obtained from the following exact solution for X, V: 

X = G+ ?l- -G+ {iuM AXr + -fAXr + MAVr) , 

T 

V = iuG-^ ?i + -G+( ^AXr - iojMAVr ) , (31) 
r 

where AX, = X(r) - X(0), AK = V{t) - V{0) . 

To obtain Z{X) = (e"'^'^) we need to average over both noise fj and over the initial 
steady state distribution of Uq. We note that the solution in ( 131]) contains U'^{t) = 
{X'^{t), V^i^r)) and this has to be expressed in terms of fj and Uq- While this can be 
done we follow a different strategy which is more convenient. 
It is easier to calculate the restricted generating function 

Z(A,f/,r|f/o) = /e-"«5(f/-f/(r))\ , (32) 

\ / Uo 

where the average is performed only over noise and for a given initial condition Uq. We 
can then obtain Z{X) = f dU J dUoPss{Uo)Z{X,U,r\Uo). The steps of the calculation 
now are as follows. We first note that, because of the (5-function constraint in (1321) .we 
can obtain Q as a quadratic function of the variables fj, U, Uq- This follows from ( fT9l) 
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l=l,N 



zu[Gt^fiL + Gtj,fin] + -F^AU 

r 



(33) 



where 



AX = X - X(0), AV = V- V{0), AU = U - U{0), and 

We replace the 5-function in fl32l) by the integral representations: — U{t)) = 
j(f^a/{2TiY^ e^-^iu-u{r)) ^Y,ere = {ai,a2, ...a2N)- We then have: 

^(A>^'-l'^o) = y^e-^-(e-(^)),,,„ 



(35) 



where E{t) = -XQ - ia^ U{t). 

We need an expression for W^ir) = (X^(r), V'^{t)) which we now obtain. We note 
that since we are using a Fourier- series representation for X(t), V(t), the correct value 
at time r is obtained from the Fourier-series by setting t = t — e, e > and taking the 
limit e — > 0. Hence we obtain: 



U 



T 



X^(t), V' (t) 



oo 

Jim J2 ( ^"^M, V^M ) e- 



For large r we note that (1/t) {u}n)e~'^'^"'^ = which follows from converting the 

summation into a integral and noting that all the poles of G^{ujn) lie in the upper half 
plane. Hence we get: 



J2 ( [G^V]"^, ^co[G^v\'^ ) e— 



n=— oo 
oo 



where 



Hence we get E{t) in ( l35l) as 

oo 



G^^, iuGfi, iujG2_i, . . . ,iiijG^ I 



E{t) = -XQ 



I e 

n=— oo 



(36) 

(37) 
(38) 

(39) 



After using the full expression for vi{u) from ( l33l) to evaluate Q in (|T9l) . we obtain 
-E(r) = So + Xl^il-^n + -^-n) where (s„ + s_„) has the following quadratic form 



2A7Z 



r 



l/nl 



(40) 
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a. 



X 



FIAU - ie 



where Fi,F2 and F3 are given by (j34|) . fl37|) and fl38l) respectively. 
Similarly, one can express sq as 
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(41) 
(42) 



Ar 



So 



(43) 



- ir^o ^o4o + "o « ^ 

2 T 

where Ao,ao, ^0 and /q are all real. 

We now first evaluate averages with respect to the Gaussian distribution given in 
f l23|) for n ^ and with distribution given in f l29|) for = 0. We get 

1 



/ (s„+s_„)\ _ 
\^ /U,Uo ~ 

Hence we get 



det(J + Ari:)A„; 
1 

^det{I + XtDAo 



exp 



:exp -ao(D ^ + AtAq) 



r 

Atl 
r 



1/0 



(44) 



(-^ ln[det(J + ArDA„)] j 

\ n=— 00 / 



X exp 



E 



ia!„(D-'+ATA„)-'a„ + ^|/„P 

2 r 



(45) 



After evaluating the required matrix inverse and determinant, we take the large r limit 
to replace all the summations over n by integrations over u. This yields 



~ e^'^^^'' exp I / — 



du 



d^dji a'^ (—uj)C (uj)a(uj) 



27r[l + TLTnX{Al3-X)r{u] 



+ A7l|/MP 



where /x(A) is given by f lTTl) and 



C(u) 



i--4X^,u'\Gl 



4X-fL(^^GtiG-j^ + 2iXujGl 



N 



(46) 



4XjLCo^Gl,Gl^ - 2iXuGIj, ^ + 4Xjnco^\GtJ' 

We see from (I41|42l) that a{uj), f{u}) are linear in AU and cr. After some algebraic 
manipulations and use of the identity ( l22l) we then arrive at the following compact 
expression for {e^)mj^^- 



where 
ifi(A) = 

with 



exp (-^a^Hi a + lAU^ H2 a + ^AU^H-^ Au) , 



II + II 



(47) 
(48) 



/i(A) 



drd 



TT 



Gi,^F2F^ + Ci,2F3Ft + Ca.iFaFj + ^2.2^3^^! 
. 1 + TlTrX{A^~X)T{uj) 
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iio A =lim— / due ^ ^ w , ^ , , ■ , 49 

and 

J?3(A) = ^(/3 + /3^), (50) 

with ux)^^-^^^±^rd. 



TT i + TlTr\{/\p-\)T{u) 

Finally to get Z{\, U, t\Uo) we substitute the expression for from (147|) into 

(!35|) and perform the Gaussian integration over a. This gives 

„r«(A) 

Z(A, f/, r|f/o) ~ , , e^^^"^^ g-i(f/-+A^-if.)ffr^(^+^fjAf/)_ ^51) 

(27r)^ A/det Hi 

Putting A = in the above expression gives the steady state distribution as 

^ ^ ^ ' ^ (27r)^ v/detifi(0) ^ ^ 

From the long time solution in (1361) it can be directly verified that 

lim {UU^) = — du dLF^Fl + dnF^Fl , 
t^oo 27r L J 

and using the fact that {UU'^) is real we see that the above equals Hi{0), consistent 
with (152]) 

Now, according ([8]), the initial and the final variables Uq and U in ( 15T]) must 

factorize, which implies [H^- H2H^^H^ - H{^H^] + [H3- HiH^^H^ - H2H{Y = 
0. Since iJi and iJs are symmetric matrices, the above condition can be expressed as 

Hs - H2H-'Hj - H~'H^ = , (53) 

Using this ( l5Ti) gives 



e^^(^) /l^ \ /l^ \ 

^(^>^--l>^°)- p,). VdetH.lA) ''"'!-^^ ^''^'^) exp(--C/„-L.(A)t/„] .(54) 

This means that we can make the following identifications for "^{U, A), x{Uo, A): 

*(?7, A) = — — —L^^^ exp (-lu^Lii\)u] , (55) 



(27r)^ v/deti?i(A) V 2 

x(f/o,A)=exp(^-i[/o^L2(A)f/o^ , (56) 

where Xi(A) = fi"f ^ + H^^H^ (57) 

and X2(A) = H2H^^HJ - = -H^^HJ . (58) 

Thus we have obtained the left and right eigenvectors of the Fokker-Planck operator Cx 
corresponding to the eigenvalue /i(A). It can be seen that the orthonormality condition 
/ dUxiU, A)^(f/, A) = 1 is satisfied. 
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We obtain Z{X) by integrating Z{X,U,t\Uq) over U and then averaging over the 
initial condition Uq with respect to the steady state distribution Pss(f^o)- This then 
gives our final expression for the correction to the CGF: 

^(A) = (^det//i(A)det/fi(0)detLi(A)det[fff^(0) + L2(A)]) (59) 
Since Li{0) = i?r^(0) and ^2(0) = 0, it is verified that 51(0) = 1. 



5. Example of Single Brownian particle 

The Langevin equation for a single Brownian particle is given by : 

mi) = -{-fL + 'yR)v + r]L{t) + r]R(t), (60) 

where v is the velocity of the particle and m is it's mass. Here we consider the velocity 
of the particle not the position since velocity v of the particle will have a normalized 
steady state distribution whereas position will not have and the heat transfer Q, in 
which we are interested, does not depend on position. For single Brownian particle the 
matrix defined in (fT3il becomes a complex number: G~^{u) = l/{—mu^ + iu'-f) where 
7 = 72, + 7r. Following all the steps described in the last section one can easily arrive 
at the expression P7|) where, //(A) is given in ffTTj) and i/'s are given in fHHjl . fH9l) and 
(150|) . In this case one can carry out the integrations present in the expressions of all 
these quantities. The expression for phonon transmission coefficient is obtained from 
(fT2l) and given by T{uj) = 47^7/2 [m^cu^ + (7L + 7r)^]~^ We use this form in (ITTi) to 
evaluate the integral and get 



l-,/l + ^:^^^T.T,A(A/?-A) 




(61) 



This is the result obtained in [16]. Similarly using G^{uj) = l/(— mw^ + iu'-f) and the 
above form for T(a;) we evaluate iJ's given by (l48l) . (149|) and (!50l) to get 

dL + dn 



Hi 



„ Adi + i(7-vW^) 



H 



a/7^ + «■ 

mA 



3 



7l + Arfi j , 



where, a = ^4:dLdR\{Af3 - A). (62) 



It is easily verified that iJs' for the single Brownian particle satisfy the relation ([53 
i.e. {H2Y + H2 — -H^i-ffa = 0. Now using the expression fl59|) we obtain 



9(A) = 2,/ ;V7^ + 4W(A/;-A) 



which also agrees with the result in 
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6. Discussions 



We have presented a formalism to calculate the CGF /i(A) and it's correction g[X) for 
heat transport across a harmonic chain connected to white-noise Langevin reservoirs. 
The formula for /i(A) is expressed as an integral over frequencies, with the integrand 
depending explicitly on the phonon transmission function T{uj). The expression for 
g{X) is in terms of integrals involving an appropriate phonon Green's function. We 
have illustrated the usefulness of the formalism by calculating yu(A) and g{X) for a single 
Brownian particle for which case all integrals can be performed explicitly. For systems 
with more number of particles, the function T{u:) can easily be obtained analytically 
for an ordered harmonic chain and numerically for disordered harmonic chains. Hence 
our formalism can be used to numerically compute /x(A) and g{X) with high accuracy. 
A knowledge of these functions would enable one to check the validity of the fluctuation 
symmetry for the large deviation function. We note that /i(A) itself is a useful quantity, 
containing information on current moments in the nonequilibrium state. We show that 
it always satisfies the fluctuation symmetry relation. Finally we have pointed out that 
yu(A) can, in general, be shown to be the largest eigenvalue of a Fokker-Planck type 
operator [Cx for our problem). Using our formalism we not only obtain this eigenvalue 
but also the corresponding left and right eigenvectors. 

The present approach has recently been generalized to the problem of computing 
/i(A) for the case of heat conduction across arbitrary harmonic networks |10]. The 
problem of calculating g{\) in such cases and also the extension of the present formalism 
to quantum systems are interesting open problems. 



Appendix A. The Fokker-Planck equation 

Let t/, t|t/o) denotes the probability distribution of heat flow Q in duration r given 
the initial and final configuration U and Uq respectively. The distribution P(Q, U,t\Uo) 
satisfies the following Fokker-Planck equation 

aP ^ r ^ d (AUi) d (AQ) 1 ^ ^ (AUiAU^) 
dt I ^ dUi At dQ At 2^^dUidU^ At 

1=1 1=1 m=l 



2N 



{AUiAQ) 1 (Ag2) 



1=1 



where the moments are calculated using the Langevin equations (jlj) and heat equation 
given in ([6]). After calculating the moments we get 

^ = CQP{Q,U,t\Uo) (A.2) 

where 

. ^ 1 [OH d m d 



i-i idxi dvi dvi dxi- 
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, IL d ^ IR d ( 2 IlTl 

H T^t-i H 7^ Vn + iLVi 



+ 



mi dvi 



nil 



_d_ 

dQ 



jrTr , . ^ ^2 , 2^lTl 



-m. (A.3) 



mf dvf dvlj- ' '^^^"''dQ'^ ' nii dvidQ 

The corresponding Fokker-Planck equation for the restricted characteristic function 
■Z'(A, t/, t|?7o) = f^^dQe~'^'^P{Q,U,t\Uo) is obtained by muhiplying both sides of the 
above equation by e~^*^ and then integrating with respect to Q. We get 

^-Z{X,U,t\Uo)=Cx Z{X,U,r\Uo) 



where 



dr 



(A.4) 



N 



on d 



1=1 



mi Idxi dvi dvi dxi 



H 7^ — vi H — Vn 



mi dvi 

^2 



H — h 



mjv dvN 
dR 



ILVI 



2^2' 2 n 2 + dlvlX^ + —^X^Vi. 

mi ovf mfj- ovf 



di 
mi 

mi dvi 



X 



(A.5) 



Appendix B. Multidimensional Gaussian integral of complex variables 

For ease of reference we give the following result for complex Gaussian integrals. 



rf"Zexp(-Z''AZ* + Z'B + C^Z*^ 



IT 



detA 



exp{C'A-'B) {B.l] 



where Z is a n-dimensional complex vector, A is a Hermitian matrix and B, C 
are arbitrary complex n-dimensional vectors. The integration J d^Z denotes the real 
integrations J d^Xd^Y with the substitution Z = X + %Y ^ X and Y being real n- 
dimensional vectors. 

Proof: Let K be the unitary matrix such that A = KDK'^ where Z) is a diagonal 
matrix with all real diagonal elements. Now with the following transformations 

Z = K'^Z , B = K^B , C = K^C 



(B.2) 



we can write the complex Gaussian integration in the form 
d'^Xd^Yexpi-Z^DZ* + Z^B + C^Z*). 



(B.3) 

Since D is diagonal, the above is a product of 2n uncoupled Gaussian integrations. 
Performing the integrations we get: 



nTT Bid . 



i=l 



detA 



exp{C' D-'B) 



detA 



exp{C'A-'B) , 



(B.4) 



which completes the proof. 
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